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Abstract 

The  sodium  borohydride  and  hydrogen 
peroxide  liquid  fuel  cell  developed  at  the 
University  of  Illinois  shows  promise  as  a  viable 
energy  source  for  a  wide  range  of  applications. 
To  achieve  higher  powers  for  a  fixed  active  area, 
an  optimal  flow  field  design  is  desired.  To  aid  in 
this  venture  a  coupled  2D-3D  model  of  the  fuel 
cell  was  developed  using  the  COMSOL 
Multiphysics  software  package.  At  this  stage  in 
development,  we  are  comparing  the  model  to 
experimental  data  to  develop  an  understanding  of 
its  predictive  capabilities.  The  model  is 
governed  by  the  Butler-Volmer  equation, 
Navier-Stokes  incompressible  flow  equations, 
Darcy’s  law,  and  the  convection  and  diffusion 
equations  that  express  the  electrical 
characteristics,  momentum  conservation,  and 
mass  conservation  respectively.  Results  show 
modifications  in  the  serpentine  flow  field  design 
can  create  a  significant  improvement  in  current 
generation.  Calculated  voltage  versus  current 
plots  for  the  standard  design  matched 
experimental  results  accurately,  but  deviations 
entered  when  other  flow  fields  were  used.  These 
inaccurates  are  mainly  attributed  to  change  in  the 
over  potential  created  by  change  in  the  flow 
field.  Despite  these  difficulties  the  model 
correctly  predicted  the  highest  performing  flow 
field  designs.  With  these  modifications  the 
power  can  be  signifigantly  increased  compared 
to  the  standard  case. 


Introduction 

The  COMSOL  model  employed  in  these 
studies  represents  the  cathode  of  a  working  cell. 
A  sketch  of  the  model  is  shown  in  fig.  1 .  The 
anode  is  introduced  through  boundary  conditions 
and  is  assumed  to  be  working  ideally.  This  is 
done  to  reduce  the  computational  requirements 
of  the  model  so  that  a  solution  can  be  found  in  a 
reasonable  time,  and  to  allow  separate 
optimization  of  the  cathode  and  anode  in  case  the 
best  configuration  is  different  for  each. 

The  first  and  last  legs  of  the  continuous 
channel  seen  in  fig.  1  are  the  locations  of  the 
inlet  and  outlet.  The  serpentine  design  is  a  one 
dimensional  representation,  but  the  assembly  is  a 
three  dimensional  structure.  While  the  reactant 
flows  along  the  channel  from  inlet  to  outlet, 
some  fluid  diffuses  into  the  diffusion  layer. 
Transport  through  the  diffusion  layer  involves 


diffusion  and  permeation  through  the  porous 
structure  of  the  layer.  The  reactant,  H2O2,  is 
consumed  and  water  is  produced  at  the  bottom 
surface  of  the  diffusion  layer.  The  current 
density  is  computed  at  this  interface. 


Figure  1:  The  Model  Domain 


Governing  Equations 

The  governing  equations  used  in  the  model 
are  taken  from  the  COMSOL  package  and 
include:  current  generation,  momentum 

conservation,  and  mass  conservation.  Each  is 
discussed  in  the  following  sections. 


Current  Generation 

Electrical  current  in  this  treatment  is 
generated  at  the  interface  between  the  membrane 
and  the  diffusion  layer,  where  the  catalyst  is 
deposited.  This  is  a  first  order  representation  of 
the  actual  cell  where  the  catalyst  is  distributed 
over  the  diffusion  layer.  Within  the  model  this 
region  exists  as  a  boundary  on  the  bottom  of  the 
three  dimensional  diffusion  layer.  The  electrical 
characteristics  are  tied  to  chemical  reaction  rates 
through  the  Tafel  type  Butler-Volmer  equation  1: 


ic  =  -(-saa„) 
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-exp 
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cH20, rf  '  ^  nRT 
Here,  the  current  density  is  a  function  of  the 
specific  surface  area  of  the  active  layer,  S„,  the 
number  of  electrons  transferred  in  the  reaction,  n. 
Other  parameters  are  Faraday’s  constant,  F,  the 
ideal  gas  constant  R,  the  cell  temperature,  T,  and 
the  overpotential,  r\op. 

While  the  reversible  potential,  which  is  the 
theoretical  half-cell  potential,  is  a  constant  for  a 
given  temperature,  the  fuel  cell  potential  varies 
due  to  voltage  losses  due  to  internal  resistance. 
In  order  to  generate  a  current  versus  voltage  plot, 
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one  of  the  main  functions  of  the  model,  the 
overpotential  must  be  known  for  a  given  cell 
voltage.  The  overpotential  function  is  estimated 
as  described  later  using  data  from  other  sources, 
including  other  models  and  experimentation. 

Momentum  Conservation 

In  the  two  dimensional  channel  system  the 
flow  is  defined  by  the  incompressible  Navier- 
Stokes  equations.  Momentum  conservation  must 
include  both  channel  flow  and  porous  flow, 
governed  by  Darcy’s  law,  in  the  diffusion  layer. 
The  two  regions  are  defined  in  separate  domains 
connected  by  a  pressure  boundary. 

i. Flow  Field  Momentum  Conservation 

The  flow  of  reactant  within  the  flow  field  is 
defined  by  the  Navier-Stokes  equations  [1,  2], 

-  V-  ^Vu +(Vu)r)+p(u .  V)u + Vpmv  =0  (2) 

V  •  u  =  0  (3) 

This  form  is  called  the  total  stress 
formulation  and  is  the  version  used  by 
COMSOL.  The  first  equation  represents 
conservation  of  momentum,  and  the  second  the 
continuity  equation.  In  this  equation,  p ,  is  the 
fluid  density,  p  is  the  viscosity,  and  u  is  the 
velocity  vector  within  the  channels.  Two 
dimensions  are  retained  in  the  channels,  whereas 
3-D  flow  is  used  in  the  diffusion  layer.  The  two 
are  linked  through  the  channel  boundary 
condition.  The  assumption  of  incompressible 
flow  is  made  in  order  to  maintain  a  constant  fluid 
density.  This  allows  the  decoupling  of  these  two 
equations  from  the  energy  equation,  assuming 
that  the  viscosity  is  constant  as  well  [3].  It  can 
be  seen  that  the  inlet-outlet  pressure  difference 
drives  the  flow  in  the  channel.  In  the  model,  an 
outlet  pressure  is  set  as  one  boundary  condition 
while  a  velocity  at  the  inlet  serves  as  the  second 
boundary  condition  to  determine  the  pressure 
difference  over  the  entire  channel.  The  inlet 
velocity  is  perpendicular  to  the  inlet  boundary. 
The  inner  walls  of  the  channel  have  a  no-slip 
boundary  condition. 

The  solution  gives  a  pressure  distribution  and 
a  velocity  profile  throughout  the  entire  channel 
system.  These  distributions  are  governed  by  the 
pressure  driven  flow  along  channels  combined 
with  flow  across  the  diffusion  layer  boundaries. 

ii.  Diffusion  Layer  Momentum  Conservation 

The  momentum  conservation  within  the 

diffusion  layer  is  defined  by  Darcy’s  law  of  flow 
through  a  porous  media  [4]. 


u 


dare 


(4) 


V  R  J 

Here  the  velocity  vector  within  the  diffusion 
layer,  u jarc,  is  proportional  to  the  permeability,  k, 
the  inverse  of  viscosity,  p,  and  the  pressure 
gradient,  Ap,  with  the  dependent  variable  being 
pressure.  Both  the  viscosity  and  permeability  are 
assumed  to  be  constant  throughout  the  entire 
diffusion  layer.  The  version  of  Darcy’s  law  used 
by  COMSOL  Multiphysics  is  shown  in  equation 


5. 


This  form  sets  the  net  flow  in  and  out  of  the 
control  volume  equal  to  the  source  term  F. 
However,  no  sources  occur  in  the  subdomain  of 
the  diffusion  layer,  only  on  the  boundaries, 
removing  the  source  term  from  the  equation. 
The  ideal  conditions  for  uniform  flow  through 
the  diffusion  layer  are  a  high  permeability  and  a 
low  fluid  viscosity,  but  this  is  not  fully  achieved 
in  practice.  There  are  two  types  of  boundary 
conditions  that  can  be  implemented  for  driving 
the  flow:  a  pressure  condition  or  a  velocity 
condition  on  a  boundary.  These  conditions  are 
both  utilized  on  various  surfaces.  As  the  Navier- 
Stokes  flow  in  the  channels  produces  a  pressure 
distribution,  these  pressures  form  a  top  boundary 
condition  for  the  Darcy’s  law  equation.  The  land 
surface  area  between  the  channels  has  an 
insulation  boundary  condition.  This  is  also  true 
for  any  borders  along  the  edges  of  the  diffusion 
layer  that  would  be  contacting  the  gasket, 
preventing  any  external  flow.  A  Flux  velocity 
condition  is  imposed  at  the  active  anode  area, 
attributed  to  the  reactant  loss  and  water  creation, 
as  well  as  water  transportation  across  the 
membrane. 


Mass  Conservation 

Mass  conservation  for  a  species  with 
concentration  c  is  governed  by  net  flow 
associated  with  diffusion  and  convection  [5]. 
The  form  used  by  COMSOL  is  called  the  non¬ 
conservative  formulation  and  is  given  in 
equation  6. 
dc 

8ls—  +  Vf-DVc)  =  Rc-uVc  (6) 
dt 

There  are  two  separate  sets  of  these 
equations,  one  set  for  H2O2  and  another  for  the 
water  forming  the  approximately  30%  peroxide 
solution.  The  time-scaling  coefficient,  8ts,  is 


zero  for  steady-state  operation,  eliminating  time 
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deviation.  Rc  is  a  reaction  rate  that  exists 
throughout  the  entire  domain.  Since  current 
generation  was  restricted  to  the  active  area 
boundary,  R,  is  also  equal  to  zero.  The  diffusion 
coefficients,  D,  are  inputted  by  the  model  user 
and  are  specific  to  the  reactant  used.  For  the 
peroxide  solution  used,  DFFO  is  le-6  m2/s  and 
DH202  is  1.37e-9  m2/s.  The  velocity  vector,  u,  is 
supplied  from  the  solution  of  Darcy’s  equation, 
thus  coupling  these  equations. 

Three  active  boundary  conditions  are 
possible.  The  first  is  a  molar  flux,  the  second  a 
concentration  and  the  third,  a  convective  flux. 
The  first  two  are  used  with  the  flux  existing  at 
the  active  area  plus  the  concentration  occurring 
at  the  boundary  shared  with  the  flow  design. 
Other  areas,  such  as  the  land  area  and  the 
diffusion  layer  edges  have  an  insulation 
boundary  condition,  at  which  no  concentration  or 
flux  exists. 

The  fluxes  defined  as  the  boundary 
conditions  are  dependent  on  the  amount  of 
electrical  current  that  is  generated  in  the  cell. 
For  every  two  moles  of  electrons  produced  at  the 
anode,  a  single  mole  of  FI202  is  consumed. 
Similarly,  a  single  mole  of  water  is  produced  for 
every  mole  of  electrons  used.  With  this  known, 
a  flux  due  to  the  reaction  can  be  found  in  relation 
to  the  current  density  by  converting  amps  per 
second  into  moles  of  electrons  per  second  by 
Faraday’s  constant. 

Water  is  also  transported  along  with  protons 
across  the  membrane  in  a  phenomenon  called 
electro-osmotic  drag.  As  protons  cross  through 
the  membrane,  they  tend  to  take  water  molecules 
with  them  into  the  cathode.  The  drag  coefficient, 
termed  drag  in  the  COMSOL  model,  is  the 
number  of  moles  of  water  entering  the  cathode 
per  mole  of  protons.  The  precise  value  of  drag 
is  dependent  on  many  factors  including 
membrane  water  content,  pressure  differentials 
between  the  anode  and  cathode,  and 
concentration  differentials  [6].  Flowever,  in  the 
present  calculations,  the  “traditional”  value  for 
transport  in  Nafion  of  drag=  3  was  assumed. 

Note  that  COMSOL  represents  molar  fluxes 
as  inward,  or  into  the  boundary.  Therefore  a  loss 
of  reactant  from  the  domain  is  actually  a  positive 
value.  Flence  the  FL02  flux  is  positive, 
representing  a  loss  of  reactant,  and  vice  versa  for 
water,  representing  a  gain. 

This  completes  the  boundary  conditions  for 
mass  conservation  and  the  flow  model.  As 
witnessed,  the  model  is  highly  coupled,  with 
information  shared  in  many  ways  across  model 
domains  via  boundaries. 


Geometry  and  Extrusion  Coupling  Variables 

The  half-cell  model  consists  of  two  separate 
geometries  that  interact  with  each  other:  fuel  cell 
flow  channels  which  are  two  dimensional,  and 
the  diffusion  layer  which  is  three  dimensional.. 
Both  are  easily  constructed  using  only  a  minimal 
number  of  geometric  functions  provided  by  the 
COMSOL  Multiphysics  solid  modeling  menus. 

By  default,  these  two  separate  geometries  do 
not  automatically  share  information  from  their 
respective  physics  modes.  Instead,  the 
information  to  be  shared  must  be  manually 
selected.  This  is  the  case  for  the  pressure 
coupling  between  the  channels  and  diffusion 
layer.  To  achieve  a  proper  setup,  the  COMSOL 
function  called  “extrusion  coupling  variables” 
was  implemented.  This  allowed  transfer  of  data 
from  a  subdomain  to  a  boundary.  For  the  present 
model,  the  expression  transferred  is  p,  and  the 
name  it  was  given  was  p_nav.  This  name  was 
then  used  in  the  boundary  conditions  to  create 
the  pressure  boundary  condition.  The  flow 
between  the  channels  and  the  diffusion  layer  is 
also  a  boundary  condition  and  is  described 
elsewhere.  With  this  in  place,  the  model  was 
complete  and  ready  to  be  solved. 

Solution  Procedure 

Initial  attempts  to  solve  the  model  as  a  single 
entity  from  the  initial  conditions  resulted  in  non¬ 
convergence.  It  was  found  that  the  couplings 
were  too  complex  for  the  modeling  software  to 
evaluate  simultaneously  without  a  fairly  accurate 
initial  value.  Therefore  a  method  to  slowly  build 
a  better  initial  value  was  developed. 

There  are  three  steps  to  the  process.  The  first 
step  is  to  remove  the  coupling  between  the  mass 
conservation  and  the  momentum  conservation. 
The  two  momentum  conservation  modes  are  then 
solved  together  using  first  rough  solution. 

With  the  momentum  conservation  solved  for, 
this  partial  solution  can  then  be  used  to  solve  for 
the  mass  balance.  By  reinserting  the  coupling, 
the  velocity  profile,  found  without  the  velocity 
due  to  the  species  flux,  is  stored  in  the 
convection  and  diffusion  mode.  Mass 
conservation  for  each  species  is  then  solved 
individually.  For  the  velocity  profile  to  exist  the 
previous  solution  must  be  used  as  the  initial 
value  in  the  solver. 

Finally,  using  the  decoupled  solutions  to  all 
of  the  physics  modes,  the  momentum 
conservation  dependence  on  the  mass  flux  is  re¬ 
entered,  and  a  solution  for  the  entire  model  using 
the  previous  solution  as  the  initial  value  is 
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carried  out.  This  process  initially  requires  a  low 
overpotential,  near  //= 0.5,  to  work  effectively. 

As  other  overpotentials  are  required  in  order  to 
develop  the  desired  voltage  versus  current  plot, 
the  overpotential  can  be  gradually  increased, 
using  the  old  solution  at  the  previous 
overpotential  as  the  initial  value.  This  also  holds 
true  for  the  variation  of  other  parameters,  such  as 
the  exchange  current  density.  Following  this 
process  provides  a  well  defined  and  smooth  V-I 
curve. 

Model  Normalization  and  Results 

With  all  parameters  evaluated  and  the  model 
constructed,  the  model  is  used  to  generate  a 
current  profile  within  the  cell.  As  an  initial  test 
to  ascertain  that  the  model  was  behaving 
properly,  the  model  was  run  at  a  base 
overpotential  ijop  =  0.5  [7]. 

While  it  was  impossible  to  determine  if  the 
average  electrical  current  density  was  correct 
until  the  result  was  fully  benchmarked  against 
experiments,  it  was  still  possible  to  judge  the 
general  validity  of  the  predicted  current 
distribution.  The  expected  distribution  shoidd 
follow  the  concentration  profile,  with  higher 
currents  directly  underneath  the  flow  field,  and 
lower  current  densities  underneath  the  land  area. 

The  model  did  simulate  this  characteristic,  with  a 
typical  distribution  shown  in  fig.  2. 

Mto:  3444.789 


Figure  2:  Current  density  distribution  for  r\op  =  0.5 


Another  observation  is  that  the  regions  of 
lowest  electrical  current  density  are  those  under 
the  channels,  where  no  cross-flow  exists.  Here, 
only  permeation  drives  the  transport  of  reactant 
to  the  active  area.  Therefore  theoretically  it  is 
advantageous  to  design  the  cell  flow  field  to 
achieve  a  maximum  active  area.  While  the 
current  NaBH4/H202  fuel  cell  flow  field  does 
extend  to  the  active  area  boundary  at  the  wall,  a 
small  artificical  gap  was  necessary  in  the  model 


due  to  meshing  constraints.  This  gap  was 
included  in  all  flow  designs  studied.  This  effect 
should  be  roughly  the  same  in  all  configurations 
and  is  thought  to  be  negligible. 

Finally,  a  relation  between  cell  voltage  and 
cell  overpotential  must  be  found  to  use  in  eq.  1. 
This  is  described  next. 

V-I  Results  Using  Overpotentials  from  Two 
Dimensional  Model 

The  first  attempt  at  finding  correct 
overpotential  values  for  cell  voltages  used  a  2-D 
model  developed  at  the  University  of  Illinois  [8]. 
The  model,  designed  for  the  NaBH4/H202  fuel 
cell,  requires  the  user  to  define  the  cell  voltage, 
using  this  information  to  calculate  the 
overpotential  according  to  equation  7,  in  which 
<pm  and  (pd  are  the  membrane  and  diffusion  layer 
potentials  measured  at  the  active  layer  boundary. 

Hop  =  </>rev  +  <t>m  -  <t>d  (7) 

Using  these  values  for  overpotential,  a  V-I 
curve  was  found  for  several  flow  field 
configurations.  The  land  width,  /,  and  the 
channel  width,  c,  were  varied  for  each  design. 
The  starting  voltage  for  each  flow  model  was  1.4 
volts  as  determined  in  the  NaBH4/H202 
characteristics.  The  final  voltage  was  either  0.4 
or  0.3V,  with  the  0.4V  models  unable  to  reach 
lower  voltages  due  to  a  build-up  of  extraneous 
model  data  causing  convergence  errors.  The  V-I 
curves  found  for  all  the  flow  field  designs  are 
summarized  in  fig.  3.  The  designs  shown  were 
all  serpentine  flow  but  with  varying  channel 
widths,  c,  and  land  widths,  1.  The  curves  exhibit 
the  classic  shape  attributed  to  V-I  curves,  starting 
with  a  large  drop,  typically  due  to  activation 
losses,  with  a  linear  decline  in  voltage  due  to 
resistance  losses. 


Figure  3:  V-I  Plots  for  Varying  Flow  Designs 
(c=channel  width,  Wand  width) 


For  the  majority  of  channel  widths  c,  a  small 
land  width  /  generally  gives  higher  V-I 
performance.  Thus  the  V-I  curves  also  indicate 
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that  the  performance  strongly  depends  on  finding 
the  optimal  land  width.  Land  width  is  clearly  a 
crucial  parameter  to  optimize  the  cell  V-I 
characteristic. 

As  noted  earlier,  the  major  objective  of  this 
work  was  to  compare  computation  with  select 
experimental  results  using  different  channel  and 
land  widths.  The  results  for  current  density 
showed  a  large  amount  of  error  when  compared 
to  experimental  values  of  0.6  to  0.7  A/cm2. 
While  some  positive  difference  was  expected 
due  to  the  ideal  nature  of  the  half-cell  model,  the 
error  was  excessive  in  terms  of  absolute  current 
values.  It  is  planned  to  correct  this  in  future  work 
by  revising  the  overpotential  correlation.  Still, 
the  comparative  results  have  provided  some 
valuable  information  relative  to  localizing  the 
region  of  optimal  performance. 


previous  results.  The  updated  V-I  curves  are 
seen  in  fig.  5. 


Figure  5:  Voltage  versus  Current  using  Updated 
Overpotentials 


V-I  Results  Using  Overpotentials  from 
Experimentation 

While  some  important  factors  in  flow  field 
design  were  observed  using  the  overpotential 
values  found  from  the  two  dimensional  model, 
the  electrical  current  and  power  density  values 
found  were  significantly  different  from 
experimental.  Therefore,  an  attempt  was  made 
to  normalize  the  model  to  past  experimental 
results  using  the  standard  “reference”  fuel  cell 
flow  field  design.  The  currents  found 

experimentally  were  compared  to  currents  found 
with  the  model.  By  determining  the  voltage  that 
the  model  assigned  for  the  experimental  current 
through  a  linear  fit  between  model  points,  the 
corresponding  overpotentials  were  found.  The 
resulting  V-I  curve  proved  an  excellent  match  to 
the  experimental  results  for  the  reference  design, 
as  seen  in  fig.  4. 


Figure  4:  Experimental  and  Normalized  Model  V-I 
Curve 


Next,  several  of  the  alternate  flow  field 
models  were  modified  to  use  these  overpotentials 
in  order  to  see  any  differences  between  the 


Observations 

To  verify  the  COMSOL  mode,  the  V-I  curves 
found  with  the  model  for  non-standard  flow 
configurations  were  compared  to  experimental 


Figure  6:  Comparison  of  nomialized  model  to 
experimental  results  for  non-standard  flow  field 
configurations 


As  expected,  the  predicted  curve  for  the 
standard  configuration  closely  matches  the 
experimental  results.  However,  the  modeled 
results  for  non-standard  configurations  show 
significant  deviation  from  the  experimental  data. 
This  indicates  that  normalizing  the  model  to 
overpotentials  found  from  the  standard 
configuration  loses  accuracy  for  significant 
changes  to  flow  field  design.  However,  as 
shown  from  fig.  6,  the  model  successfully 
predicted  the  two  highest  V-I  configurations. 

Conclusion 
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An  important  aspect  of  the  computational 
technique  is  the  three  step  procedure  described 
that  obtains  an  appropriate  initial  condition  to 
allow  convergence  on  a  desktop  PC.  Another 
key  point  involves  the  overpotential  correlation 
which,  while  a  step  in  the  right  direction,  still 
needs  improvement.  This  model  is  still  a  work 
in  progress  and  does  not  yet  fully  predict  the 
performance  of  non-standard  flow  field  designs. 
However,  it  is  already  a  useful  tool  for 
determining  trends  and  the  region  of  optimal 
performance. 

It  is  vastly  preferable  to  use  modeling  to 
localize  the  region  of  interest  rather  than 
machine  new  flow  plate  designs  for  each  data 
point. 

The  most  important  design  result  found  from 
the  modeling  process  is  that  the  V-I 
characteristic  of  the  present  NaBH4/H202  cell 
could  be  significantly  improved  by  redesigning 
the  flow  field  to  reduce  both  land  and  channel 
width.  This  change  in  flow  plate  design  requires 
the  replacement  of  the  present  graphite  plates 
with  pieces  having  channel/land  dimensions  in 
the  region  of  the  optimal  space  of  fig.  6.  No 
other  changes  in  the  construction  process  or  the 
materials  used  are  required. 

In  conclusion,  while  the  model  still  contains 
inaccuracies  the  desired  optimal  flow  design  can 
be  finalized  with  a  few  flow-limited  experiments 
located  in  the  region  of  optimization  predicted. 
Further  information  can  be  found  in  [9]. 
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Appendix:  Constants 


Sa 

Specific  surface  area  of  active  layer 

108  nr/m3 

d 

Active  layer  thickness 

10"6  m 

n 

Number  of  electrons  transferred  in  reaction 

2 

F 

Faraday’s  constant 

96485  J/mol  K 

R 

Ideal  gas  constant 

8.314 

T 

Cell  temperature 

303.0  K 

^lop 

Cell  overpotential 

0.5  V 

k 

Permeability 

M011  nr 

cH202 

Concentration  of  F1202 

3.029e+3  M 

dh2o2 

Diffusion  coefficient  of  H202 

1.366e-9  m2/s 

cH20 

Concentration  of  F120 

5.150e+4M 

DFEO 

Diffusion  coefficient  of  H20 

le-6  m2/s 

Rc 

Reaction  rate 

0 

Drag 

Electro-osmotic  drag  coefficient 

3.0 

P 

Fluid  viscosity 

1.056- 10”3  Pa-s 

P 

Fluid  density 

1003  kg/m3 
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Aft- 


UIUC 


Equipment  for  Testing  of  500W 
Design  includes: 


N1  DAQ  (Labview) 
Programmable  Load 
2  KW  Power  Supply 
Water  Deionization 
Hot  Press 

Pumps  /  Electronics 


Screen  Printer 
Digital  Oscilloscope 
Hydrogen  Detectors 
Fume  Hood 
Safety  Shower 


Fuel  Cell  Testing  Lab 


6/22/2006-6/24/2006 


COMSOL  Users  Conf.  2006  Boston, 
MA 


3 


NaBH4/H  O  Fuel  Cell 


♦  Use  in  fuel  cells  is  a  relatively  new  development 

♦  H2/H202  and  NaBH4/H202  cells  were  investigated  at 
NPL  Associates,  Inc.,  the  University  of  Illinois 
(UIUC),  and  elsewhere 

♦  Have  shown  great  results,  demonstrating  the 
general  feasibility  of  a  peroxide  based  cell 

♦  Excellent  potential  for  space  applications  due  to 
high  power  density  and  air  (oxygen)  independence. 
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UIUC/NPL  Direct  Peroxide  Fuel  Cells 


♦  The  sodium  borohydride/hydrogen 
peroxide  reactions. 

Anode: 

NaBH4  +  4 H20  ->  NaB02  •  2 H20  +  8 H+  +  8e” 

Cathode: 

H202  +  2 H+  +  2e  ->  2 H20 
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Peroxide-based  fuel  cells  (both  the  H2/H202  and  NaBH4/H202  types) 
offer  higher  V-I  densities  than  traditional  fuel  cells  under  similar 
conditions 


Performance  Comparison  of  Various  Fuel  Cells 


The  V-I  characteristics  of  various  fuel  cells,  at  room  temperature,  ambient  pressure  operation. 

COMSOL  Users  Conf.  2006  Boston, 

MA 
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II  -  MEA  studies  at  UIUC 


♦  The  UIUC  Reference  MEA 

♦  Tests  with  Various  Catalysts  (anode  and  cathode) 

♦  Results  from  the  10-40  W-class  Test  Fuel  Cells 

♦  Implication  for  the  500-W  Stack 

♦  Initial  500-W  Stack  testing. 
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10-40  Watt  Class  Catalysis  Test  Cell 


■  5  x  5cm  flow  field/active  area 

■  Nafion  115  untreated  w/o  catalyst 

■  Nafion  binder  for  catalyst  on 
diffusion  layer 

■  Graphite  felt  diffusion  layer 

■  Graphite  current  collectors 

■  Aluminum  endplates 

■  Silver  epoxy  at  interface  between 
endplates  and  graphite  plates 


10  Watt  class  fuel  cell 
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Test  Cells  -  Compact  10-40  W  Power  Units 


The  15-W  cell  shown  here  uses  an 
integrated  cooling  channel  to  dissipate 
the  waste  heat  generated  in  the 
relative  small  25-cm2  active  area. 

An  optimized  version  of  this  small  cell 
generated  36-W  at  ~  60°C, 
representing  the  highest  power 
density  reported  to  date  for  a  small 
fuel  cell  working  at  sub-1 00°C. 


15-W  NaBH4/H202  Test  Fuel  Cell  as  assembled. 
■  Flow  rate  of  approximately  200  cm3/min 

■Minimal  pressure  drop  even  with  parallel  flow  due  to  low  flow  rate 

■Temperature  rise  of  approximately  15°C 

■Heat  flux  is  approximately  equal  to  electrical  power  (500-W) 
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Unique  Catalyst  Issues 


♦  Key  issue 

■  Our  reactants  decompose  violently  (gas  production)  in 
the  presence  of  conventional  fuel  cell  catalysts  such  as  Pt 
as  well  as  several  other  very  active  catalysts 

■  This  dramatically  reduce  the  cell  efficiency 

♦  Goals  of  catalysis  studies 

■  Find  catalysts  for  both  anode  and  cathode  reactions  that 
maximize  power  density  without  gas  production 
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Voltage  (V) 


UIUC/NPL  Proprietary  Catalysts 


0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6 


Catalysts  optimized  for  high  efficiency 
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Power  (W) 


500-W  Stack  Catalyst  Selection 


♦For  initial  tests  of  the  500-W 
stack,  palladium  and  gold 
alloys  were  finally  selected  for 
the  anode  and  cathode 
catalysts,  respectively,  for 
their  high  performance  and 
relative  stability 
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The  500-W  UIUC/NPL  NaBH4/H202  Fuel  Cell  Stack 


The  active  area  per  cell  was  144  cm2  and  15  cells  were  employed  to  provide  a  total  stack  active  area  of 
2160  cm2. 
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Bipolar  Plate  Design 


□  Serpentine  design  utilized 

□  Ensures  proper  fuel  distribution  over  range  of  environments 

□  Higher  pressures  caused  by  serpentine  are  within  design 
parameters 

□  External  Inlets  and  Outlets 

□  NaBH4  conductivity  causes  risk  of  short  circuit  when 
manifolding 

□  Non-conductive  tubing  between  cells  during  series  flow 
prevents  short  circuit 
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Stack  Parameters 


□  Flow  rate  ~  200  cm3/min 

□  Pressure  drop  for  series  flow  is  17  psi  for  cathode  and  8 
psi  for  anode 

□  Expected  temperature  rise  ~  15C 

□  Heat  flux  ~  500  W 
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500-W  Stack  runs  at  ~  29  A,  17  V 


Performance  of  NPL  500W  NaBH4/H202  Fuel  Cell 


Current  (A) 


V-l  and  P-l  performance  of  the  500W  UIUC/NPL/CUA  NaBH4/H202  FC. 
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Power  (W) 


Comments  re  applications:  20-W I- 
Charger™  Laptop  Power  Unit 


NaBH4/H202  I- 
Charger™  Micro  Cell 

.  20  W 

■  500  W-hr 

■  ~2.2  lbs 

SOA  Micro  DMFC 

■  20  W 

■  500  W-hr 

■  3.5 -7  lbs 
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1  -kW  I-Charger™  UPS  System 

Initial  unit  designed 
using  retail 
components. 

Compact  Design 
providing  ~1  hr 
runtime. 

Easily  expandable 
with  auxiliary  tanks 
for  extended 
operation. 

Meets  or  exceeds 
typical  commercial 
UPS  system 
requirements. 

-25%  more  space 
efficient  than  typical 
UPS 
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Objectives  for  COMSOL  modeling 


♦  Gain  insight  into  behaviors  governing  flow 
and  current  distributions 

♦  Determine  space  (diffusion  layer  parameters, 
conductivity  effects,  flow  channel  and  land 
dimensions)  for  detailed  optimization  physics 

♦  Guide  future  design  improvements 
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Model  developed  from  COMSOL 

built-in  library 


♦  Adapted  3D  Serpentine  Flow  model  for  a  cathode 
in  the  COMSOL  Model  Library. 

■  The  model  inputs,  such  as  exchange  current  density, 
concentration,  and  diffusion  coefficients  were  altered. 

■  Flow  in  the  channels  was  condensed  from  3D  to  2D. 

■  Some  physics  modes  were  replaced  with  alternate  modes 
that  best  represented  the  cell. 
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Model  Development  and 

Operation 


♦  End  result  has  a  2D  flow  field  acting  as  a  boundary 
on  the  3D  diffusion  layer. 


MA 
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COMSOL  functions  used  in 
NaBH4/H202  fuel  cell  model 


♦  Physics  Modes 

■  Momentum  Conservation 

•  Navier-Stokes  Incompressible  Flow  within  the  flow 
field 

•  Darcy’s  Law  for  flow  in  a  porous  media  within  the 
diffusion  layer. 

■  Conservation  of  Mass 

•  Chemical  Engineering  Module  Convection  and 
Diffusion  within  the  diffusion  layer. 
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COMSOL  functions  used  in 
NaBH4/H202  fuel  cell  model 


♦  Current  Density  and  Mass  Flux  Simulation 

■  Tafel  type  Butler- Volmer  Equation  used  to  find 
current  density. 

•  Depends  on  concentration  distribution  in  the  diffusion 
layer  and  over-potential  (directly  related  to  cell 
voltage). 

■  Mass  flux  at  active  area  proportional  to  the 
current  generated,  with  H202  consumed  and  H20 
generated 
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COMSOL  Application  Mode 

Coupling 
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Relationship  between  COMSOL 
functions  in  the  model 


♦  Couplings 

■  Pressure  coupling  between  Navier-Stokes  in 
flow  field  and  Darcy’s  Law  in  diffusion  layer 

■  Mass  flux  at  the  active  layer  causes  small 
velocity  boundary  condition  for  Darcy’s  Law 
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Methodology  -  convergence  requires 
special  method  to  obtain  reason  initial 

values  


♦  Solving 

■  Initial  values  found  by  first  decoupling  the 
momentum  and  mass  conservation  and  solving 
separately. 

■  The  initial  value  was  used  to  find  the  solution 
with  couplings  reintroduced 

■  Over-potential  values  were  varied  to  generate  V- 
I  curves. 
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Methodology 


♦  Solutions  were  found  for  several  different 
flow  fields. 

■  Flow  field  designs  varied  in  channel  and  land 
width. 

■  All  couplings  were  verified  and  all  boundary 
conditions  found  to  be  satisfied. 
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Initial  Results 


♦  Current  generation 
greatest  under 
channels  where 
concentration  greatest. 

♦  Current  generation 
underneath  the  land 
follows  the  flow  of 
reactant  created  by 
pressure  differences 
between  channels 

Current  Density  Distribution 
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Apparent  optimal  flow  field  dimensions 
found  by  varying  channel  and  land  widths 


♦For  most  channel  widths,  small  land  width  is 
desired. 

♦  V-I  curves  also  indicate  that  performance  greatly 
depends  on  finding  the  optimal  land  width. 
Without  a  properly  sized  land  width,  the  fuel  cell 
will  perform  well  below  potential. 
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Smaller  land  width  results  in 
greater  current  density 


V-l  Curve  for  various  flow  field  dimensions  (c=channel  width,  l=land  width 
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Model  shows  variations  in  pressure  not 
as  significant  as  land  width 


♦  By  plotting  the  current  generated  at  0.5  V  as  a 
function  of  the  overall  pressure  drop,  it  was 
determined  that  there  is  no  strong 
dependence  of  current  density  on  pressure 
drop. 

♦  Even  when  the  curve  is  modified  to  take  into 
account  variation  in  land  width,  only  a  weak 
dependence  on  pressure  is  shown. 
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Current  density  found  to  be  not 
directly  related  to  pressure 


Average  Current  Density  vs.  Model  Pressure  Drop 
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Normalized  model  V-I  curve  closely 
resembles  real  data  for  standard  cell 

♦  The  model  was 
normalized  to 
experimentally  found 
over-potentials  from 
the  standard  cell 
configuration 

♦  The  resulting 
normalized  curve  fit 
well  against  the 
experimental  data 


Experimental  and  normalized 
model  V-I  curve 
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Normalized  model  used  to  make  new  V-I 
curves  for  alternate  flow  fields 


♦  Normalized  model  used  to  generate  V-I 
curves  for  a  variety  of  alternate  flow  field 
designs 

♦  Use  of  the  normalized  model  should  result  in 
better  prediction  of  experimental  results 
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Model  indicates  1.5mm  channels, 
1mm  lands  should  perform  best 


Voltage  versus  Current  using  Updated  Overpotentials 
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Verification  of  model  with 

experiment 


♦  6  new  flow  field  plates  were  constructed 
matching  the  configurations  explored  in  the 
model 

♦  Cell  performance  was  measured  with  the  new 
flow  plates  and  compared  to  the  model’s 
predictions 
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■  Normalized  model  does  not  match  experiment 
well,  but  predicts  optimal  performance  region  in 

agreement  with  experiments. 


Comparison  of  normalized  model  to  actual  results 

for  alternate  flow  fields  -  note  that  c= 1.5mm, 

1=1  mm  performs  best,  as  predicted  by  the  model 
COMSOL  Users  Conf.  2006  Boston, 

MA 
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Normalized  model  predicts  best  design  region-  i.e. 

key  trends  predicted  are  valid 


♦  The  model  results  did  not  accurately  match 
cell  performance  away  from  the  standard 
design,  but  it  successfully  predicted  the 
design  region  for  best  configuration  shown  in 
experiments. 
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Further  lessons  learned  from  the 

model 


♦  Unlike  the  model,  experimental  results 
indicate  that  current  density  is  dependent  on 
an  ideal  combination  of  both  channel  and 
land  width,  not  just  land  width 

♦  Normalization  from  the  standard  cell  did  not 
apply  to  V-I  curves  from  alternate  flow  field 
designs,  indicating  that  the  overpotentials 
change  with  the  flow  field 
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Conclusions 


♦  COMSOL  can  be  used  to  model  the  behavior  of  an 
NaBH4/H202  fuel  cell 

♦  Model  over-potentials  must  be  normalized  to  experimental 
values  to  predict  cell  performance 

♦  Normalized  model  correctly  predicted  which  flow  field 
designs  would  work  best,  but  did  not  accurately  model 
their  performance 

♦  Normalization  is  dependent  on  flow  field  parameters,  so  a 
more  accurate  normalization  method  is  desired  in  the 
future. 
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Conclusions  -  continued 


♦  The  model  is  requires  further  work  to  accuracy 
predict  absolute  performance  values  for  case  away 
from  the  “standard”  case  used  in  normalization. 

♦  Still  it  provides  a  very  valuable  tool  for  identifying 
the  design  space  where  optimal  performance  is 
expected.  Thus  experiments  can  be  rapidly  focused 
on  this  space  to  refine  the  optimal  design. 
Significant  improvements  have  been  achieved  in 
this  fashion. 
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Thank  Y  ou 


For  more  information  please  contact: 

Dr.  George.  H.  Miley 

UIUC 

Phone:  (217)  333-3772 
email:  ghmiley@uiuc.edu 
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